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Abstract 



' We consider a one-dimensional "gas" of inelastically colliding particles where ki- 

netic energy is dissipated by the excitation of vibrational degrees of freedom. In our 
model the coefficient of restitution is a stochastic quantity whose distribution can be 
, calculated from an exact stochastic equation of motion. We investigate the equiparti- 

' tion properties of the system and propose a new algorithm for computer simulations, 

, that is a combination of event-driven(ED) and Monte-Carlo methods. 



Numerical and theoretical approaches to the dynamics of granular materials frequently 
adopt the concept of a coefficient of restitution that determines the energy loss during colli- 
sions of granular particles. Event-driven (ED) simulations ||l], |2|, ^] have shown that model 
I systems with fixed coefficient of restitution evolve into clustered states where a hydro- 

dynamic description ceases to be correct: Fundamental assumptions of hydrodynamics 
O ■ concerning the validity of molecular chaos and local equilibrium are violated Q]. On the 

other hand, molecular-dynamics simulations Q have the difficulty that ad hoc assump- 
tions about microscopic interaction laws have to be made. An inadequate choice of the 



■ interaction parameters can lead||y] to spurious effects in the simulations. 



Here we present a one-dimensional model where colliding particles interact via an ex- 
ponential potential. A control parameter is introduced that allows to perform the limit 
of a hard core interaction. On collision particles loose kinetic energy by the excitation 
of internal vibrational modes. In a statistical description we characterize the energy con- 
tained in this bath of internal oscillators by a bath temperature T^ath- Our main aim is 
to investigate the cooling properties of the system: Starting from an initial distribution of 
center of mass velocities we would like to understand in detail, how energy is transferred 
from the translational to the internal degrees of freedom. 

The principal results of this Letter are as follows: (1) Granular cooling can be un- 
derstood as energy equipartition among all (translational and vibrational) modes of the 
system. (2) At a single collision the coefficient of restitution e is a stochastic quantity. We 
derive an exact stochastic equation of motion for the relative coordinates, from which e 
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can be calculated. For the special case of "cold" particles (Ti,ath =0) we recover a clas- 
sical result of the wave theory of impact. (3) We give a statistical interpretation of the 
probability density of e and propose a new algorithm for simulations of granular particles. 

The model 

We consider a one-dimensional system of N elastic rods with mass density p and elastic 

module E. The elastic state of the ith particle at time t is described by the longitudinal 

strain field u{s,t), — ^ < s < ^, where s is the internal position coordinate and li is the 

length of the one-dimensional rod. In order to seperate internal and translational degrees 

ii 

of freedom the strain field contains only fluctuations, i.e. /^;. ds Ui{s,t) = 0, while 

homogeneous distortions contribute to the center of mass coordinate Ri of the particle. 
On collision adjacent particles interact via a potential Vhcir), which depends only on 
the current end-to-end distance of the distorted rods r = Ri^i^i + nj+i(— t) — t), 
where Ri+i^i := Ri+i— Ri— ''^~^2~^^ . In the present calculations we investigate an exponential 
interaction potential with characteristic length scale ^ 

VUr) = Be-^' , (1) 

which includes for large a the limit of a hard core potential. The constant B can be 
chosen arbitrarily, since changing its value merely corresponds to a rescaling of the relativ 
coordinates We expand the strain fields Uj in sine and cosine normal modes with 

amplitudes q^'^\t) and normal frequencies u;i^ = i'-K-^, where c=^l^ is the speed of sound. 

If we denote conjugate moments by Pi and p\'^\ respectively, and treat the elastic energy 
in harmonic approximation, the Hamiltonian reads 0: 

'H = Hbath [pf^ , ql"^ } + Ti-trans + Ti-int 

1 = 1 u=l \^ ' J ?,= 1 ' 

2 = 1 I \ V 

Here we assume periodic boundary conditions, i.e. place the particles on a ring with 
perimeter P, so that R^+i^n ■= P + Ri — Rn and q^^\^i ■= q^i^ ■ Because of the nonlin- 
ear coupling between oscillator modes our model has some resemblance with the classical 
Fermi-Past a-Ulam (FPU) problem |^]. The important difference is that here the inter- 
action is "switched" on and off depending on the relative distance of the translational 
coordinates: As long as a particle is not involved in a collision, all its vibrational modes 
are effectively decoupled. 

Equipartition properties 

We now confine ourselves to two-particle collisions and consider the case N = 2: In 
the center of mass frame we introduce the reduced mass p and an effective length scale 
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/ = 2/1^2/(^1 + ^2)- Taking N^odes normal modes for each particle we numerically integrate 
the equations of motion resulting from the Hamiltonian (|2|). The following scenario is 
considered: We start with an initial relative velocity i22,i(0) = — Vq and place the particles 
at maximum distance, i?2,i(0) = P/2 ^ ^. We concentrate on the hard core limit a ^ y^l 

and set the initial temperature of the bath of oscillators equal to zero: qf\^) =pf\Q) = d. 
The system now undergoes an infinite series of collisions with the two particles colliding 
alternately at either ends. We observe a rapid decrease of Etr = ^R2i-, the kinetic energy 
of the relative coordinate. Fig. 1 displays a decay of Etr over a range of approximately 
100 collisions until a stationary state is reached with the kinetic energy fluctuating around 
an equilibrium value /2{2N modes + 1); that is, of course, essentially zero for large 
^modes ■ We intcrprete this phenomenon as a spreading of the translational energy among 
all degrees of freedom of the system. 

This view is supported by a statistical analysis of the distribution of vibrational en- 
ergies: Denoting modal energies by E^^'^^ we define the energy of the oscillator bath as 
Ebath = 2J2u ^i'^^ ■ After a transient time we calculate the relative modal energies 
wl'^^ = e\'^'^ / Ei,ath after each collision and take the average. (Due to the continuous inter- 
change of energy between transational and vibrational modes Ei,ath itself is a dynamical 
quantity.) The resulting energy spectrum (inset of Fig. 1) shows a homogeneous partition 
of vibrational energies. This is what we shall call modal equipartition regime in the follow- 
ing. (For particular values of 7 = /i 7/2, the length ratio of the rods, we found a periodic 
structure underlying the spectra. This subtlety does not affect the main results of this 
section, however.) We have also investigated the distribution function of w^'^^ {i and v 
fixed) and always observed an excellent agreement with the Boltzmann distribution. 

An important question is, how fast the modal equipartition regime is reached. Does 
the relaxation time (measured by the number of collisions) scale with Nmodes'^ A sen- 
sitive probe of equipartition is the normalized spectral entropy r] = 1 — h/hmax, where 

h = — Yl,i V '^Y^ InzUj-'^^ and hmax = ln(2A'^modes)) which was first introduced in numerical 
studies on the FPU problemQ. If equipartition takes place, ry will decay from its ini- 
tial value T/ = 1 to its equilibrium value r/eg (for Boltzmann distributed wf^ one finds [ 10 1 
ryeg ~ 0.423/ ln(2A'^modes))- Our plot of 7? for different numbers of modes (Fig. 2) suggests a 
relaxation time of approximately 20 collisions, essentially independent of N^odes- Compar- 
ison with Fig. 1 shows that the time scale for modal equipartition is considerably shorter 
than the relaxation time of the granular cooling process. 



Stochastic coefficient of restitution 

We are now going to discuss in detail a single two-particle collision in the modal equipar- 
tition regime. Let us denote the precollisional and post-collisional relative velocity of the 
particles by —V and V' , respectively. In our model the coefficient of restitution e = V'/V 
is a stochastic quantity (cf. Fig. 1) and the collision process is described by a stochastic 
differential equation. If modal equipartition is valid, we can apply the usual concepts of 
statistical physics and consider the precollisional values of the internal coordinates as inde- 
pendent, Boltzmann distributed variables (with bath temperature Thath = Ebathf^^^modes)- 
Note that this choice of initial conditions is equivalent to imposing a Markov approxi- 
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mation on the collisional dynamics, since correlations with preceeding cohision events are 
neglected. 

Obviously, the coefficient of restitution contains the complete information about the 
energy transfer between Etr = ^V"^ and Ti^ath'- Due to conservation of total energy we have 
E'-f.j. = e^Efr and T^'^^/j = Ti,ath + 2N~'^a after collision. In order to derive the equation of 
motion for the translational coordinate it is convenient to introduce the following scaled 
variables: t = jt (scaled physical time), z(t) = ai?2,i5 = z{tq) (initial distance), and 
K = al^ (=— ^(tq), scaled initial velocity). Again we are mainly interested in the hard 
core limit so that now k;^>1 is our large control parameter. Furthermore we set the free 

2 

constant B = fij^j^. The internal coordinates in Hamilton's equations of motion can be 
integrated out0 using Green's function method and we arrive at the following equation 

dz 
dr 



for W{t) = ^ + k, the velocity increase throughout collision: 



'^^ exp {kt-zo- W{t) - ^ ^ T^(r - nF,) + Qi(r) + Q2{t) \ , (3) 

1=1,2 n=l 



dr 



with W{t) := for r < tq. The coefficient of restitution is then obtained from the 
asymptotics of W: 

e = lira - 1 . (4) 

In Eq. (^) Ti = 1 + 7 and r2 = 1 + I/7 are the periods of the lowest normal modes in 
scaled time and Qi,Q2 are Gaussian stochastic processes resulting from the initial condi- 
tions. They take on the form of stochastic Fourier series Qi(r) = J2rv=i ^ sin(n27rT/Fj) + 
J2'^=i ^ cos(n27rr/Fj) with independent and identically distributed Gaussian coefficients 
ttinjbin ■ As a consequence, we have 

< Q,{t) > = (5) 

Thnih F,; 



E 
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for < T — r' < Fj and periodic continuation everywhere else. The periodicity of the 
correlation function and the memory terms in Eq. (^) reflect the fact, that we consider 
a system of ideal oscillators that never forgets its collision history. Note that typical 
fluctuations of the noise processes are of the same order as W{t)^ which makes an analytical 
treatment even more difficult. 

However, it is possible to solve the deterministic equation (Tf,ath = 0): In the large 
K regime the coefficient of restitution is essentially given by the length ratio 7: e = 
min(7, 1/7) + ln(K|F2 — Fi|)/k, and the actual collision takes place in the time intervall 

^ +min(Fi, F2)], i.e. the duration of the collision is equal to the wave propagation 
time of the shorter rod. These results are in agreement with the wave theory of impact , 
that analyzes the collision problem in terms of stress waves propagating through the rods. 
To our knowledge, a derivation in the framework of classical mechanics (with correction 
terms for finite k) has not been given before. 

For Tbath 7^0 we had to integrate Eq. (^ numerically and average over many realizations 
of the noise processes. We have focussed on the statistics of e^, since it is the relevant 
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quantity for energy transfer. In the hard core regime (k larger than ~30) the probability 
density of e^, p(e^), does not significantly depend on k, but it is still a function of Etr/Ti^ath 
and 7. Recall that Eq. (^) describes a single collision, given the current values of Etr and 
Tbath- After collision the energy ratio has changed to £'|^/T^'^^^ (see above) and in turn 
there is a different distribution of at the next collision. Fig. 3 displays plots of the 
distribution function for different values of the parameters. 

The important point is that p{e^) has to be interpreted as a transition density: In the 
context of the Markov approximation leading to Eqs.(^),(^), the time evolution of Etr (cf. 
Fig. 1) is simply a Markovian jump process in discrete "time" (i.e. cumulative number 
of collisions). If the system is in a state with energy Etr, the probability density for a 
transition to energy E[^ is determined by 

p{Etr^El) = ^p(e^ = ^\ . (6) 



Etr \ Etr / 

As expected, the transition density resulting from Eq.(^ can be shown numerically to 
yield a stationary energy density p'^^°'^{Etr) oc exp{—Etr/Tf,ath) ■ In view of the many- 
particle problem (see below) it seems reasonable to replace the exact distribution of 
by a simpler one, that obeys detailed balance and leads to the same stationary density 
p^^"'^{Etr)- A possible choice, for example, is the analog of Glauber dynamics [ 12 1 for Ising 
spin systems: 

(G)f 2\ ^tr ( 2 ^tr \ /„x 

p^^^(e) = - exp -e . (7) 

hath \ ^ bath / 



Many— particle problem 

For the many-particle system translational energy and bath temperature become func- 
tions of time and position. From the preceeding sections it is clear that each two-particle 
collision tends to establish a state where locally Etr ~ T^ath- As in the case N = 2, we 
assume modal equipartition for each particle and treat collisions in the Markov approxima- 
tion. Furthermore, we expect that the cooling properties of the system will not sensitively 
depend on the detailed form of the probability density p(e^). 

For computer simulations we therefore suggest an algorithm that is a combination of 
ED and Monte Carlo methods: Each of the particles is assigned an individual bath 

(i) 

temperature T^^^^ that characterizes its vibrational state. For each collision the following 
steps have to be performed: (1) Calculate the next collision event by the ED procedure. 
(2) For the colliding pair determine Etr and an averaged local bath temperature T^ath- (3) 
With these parameters choose e randomly from a simple distribution, e.g. p^'^\e^). (4) 
Update the relative velocities and the post-collisional values of T^ath'-i §° back to (1). 

Another possible extension of our study would be a generalization of the inelastic 
Boltzmann equation [Q, |l| to a stochastic coefficient of restitution. Work along these lines 
is in progress. 
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Figure captions 



Figure 1: Time evolution of the relative kinetic energy for a two-particle 
system with length ratio 7 = I1/I2 = 0.6173, aVo/cZ = 4000, and N^odes = 75. 
The dashed line is the equilibrium value resulting from energy equipartition. 
The inset displays a cutout from the corresponding spectrum of modal energies: 
After a transient time of 500 collisions the relative energies wf"^ have been 
averaged over another 4000 collisions. Grey and black bars correspond to 
particle 1 and 2, respectively. 



Figure 2: Time behaviour of the spectral entropy rj for different numbers of 
modes. Again, 7 = 0.6173 and qVo/c/ = 4000. 



Figure 3: The probability density p(e^) (dotted line) and its integral, the 
distribution function (full line), for different energy and length ratios: a) 
Etr/Tbath = W, 7 = 0.6173, b) Etr/Tbath = l, 7 = 0.6173, and c) Etr/Tbath = l, 
7 = 0.2173. Data were obtained from the numerical integration of 8000 real- 
izations of Eq. (Ill) in the hard core regime («; = 100). 
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